#include <iostream>
#include <bitset>
#include <vector>
#include <fstream>
#include <string>
#include <cmath>
#include <cstdlib>
#include <cfloat>
#include <bits/stdc++.h>
#include "matrix.h"
#include "mpi.h"


using namespace std;

#define M_string "aes_"
#define FUNC 6
#define epsilon 0.001

double costfun_min[13] = {0, SIZE + epsilon, 2* SIZE + epsilon, epsilon, epsilon, SIZE + epsilon, 2* SIZE + epsilon, epsilon, 2* SIZE + epsilon, 4* SIZE + epsilon, epsilon, 2* SIZE + epsilon};
string RANK;

int minm = 100;
int sz_minm;




class Mat{
public:
    bitset<SIZE> m[SIZE];
    Mat(){};
    Mat(bitset<SIZE> rhs[SIZE]){
        for(int i = 0; i < SIZE; i++)
            m[i] = rhs[i];
    }
    Mat(const Mat& rhs){
        for(int i = 0; i < SIZE; i++)
            m[i] = rhs.m[i];
    }
    Mat(const vector<bitset<SIZE>>& str){
        for(int i = 0; i < SIZE; i++){
            m[i] = str[i];
        }
    }
    Mat inv(){
        bitset<SIZE> ret[SIZE];
        bitset<SIZE> original[SIZE];
        bitset<SIZE> tmp;
        for(int i = 0; i < SIZE; i++){
            ret[i][i] = 1;
            original[i] = m[i];
        }
        for(int i = 0; i < SIZE; i++){
            if(original[i][i] == 0){
                int k;
                for(int j = i + 1; j < SIZE; j++)
                    if(original[j][i] == 1){
                        k = j;
                        break;
                    }
                swap(original[i], original[k]);
                swap(ret[i], ret[k]);
            }
            for(int j = i + 1; j < SIZE; j++){
                if(original[j][i] == 1){
                    original[j] ^= original[i];
                    ret[j] ^= ret[i];
                }
            }
        }
        for(int i = SIZE - 1; i > 0; i--){
            for(int j = i - 1; j >= 0; j--){
                if(original[j][i] == 1){
                    original[j] ^= original[i];
                    ret[j] ^= ret[i];
                }
            }
        }
        Mat res(ret);
        return res;
    }


    Mat trans(){
        bitset<SIZE> ret[SIZE];
        for(int i = 0; i < SIZE; i++){
            for(int j = 0; j < SIZE; j++){
                ret[i][j] = this->m[j][i];
            }            
        }
        Mat res(ret);
        return res;
    }
    
    friend Mat mul(Mat a, Mat b){
        Mat ret;
        for(int i = 0; i < SIZE; i++)
            for(int j = 0; j < SIZE; j++)
                for(int k = 0; k < SIZE; k++){
                    if(a.m[i][k] && b.m[k][j])
                        ret.m[i].flip(j);
                }
        return ret;
    }


    int mxhw_row(){
      int mx = 0;
      for(int i = 0; i < SIZE; i++)
        mx = mx < m[i].count()? m[i].count(): mx;
      return mx;
    }

    int mxhw_col(){
      int mx = 0;
      for(int i = 0; i < SIZE; i++){
        int cnt = 0;
        for(int j = 0; j < SIZE; j++)
          if(m[j][i] == 1)
            cnt++;
        mx = mx < cnt? cnt: mx;
      }
        
      return mx;
    }
    
    Mat row_i2j(int i, int j){
        Mat ret = *this;
        ret.m[j] ^= ret.m[i];
        return ret;
    }
    
    Mat col_i2j(int i, int j){
        Mat ret = *this;
        for(int k = 0; k < SIZE; k++){
            if(ret.m[k][i])
                ret.m[k].flip(j);
        }
        return ret;
    }
    
    
    double cst_func(int i){
        double ret;
        switch(i){
            case 1: ret = this->h_sum(); break;
            case 2: ret = this->H_sum(); break;
            case 3: ret = this->h_prod(); break;
            case 4: ret = this->H_prod(); break;
            case 5: ret = this->h_sq(); break;
            case 6: ret = this->H_sq(); break;
            case 7: ret = this->h_prod() + this->h_prod_t(); break;
            case 8: ret = this->h_sq() + this->h_sq_t(); break;
            case 9: ret = this->HT_sq(); break;
            case 10: ret = this->HT_prod(); break;
            case 11: ret = this->H_my(); break;
            case 12: ret = this->H_var();break;
        }
        return ret;
    }
    
    void print(){
        for(int i = 0; i < SIZE; i++){
            for(int j = 0; j < SIZE; j++)
                cout << m[i][j] << ' ';
            cout << endl;
        }
        cout << endl;
    }
    
    double h_var(){
      double ret = 0;
      int a[SIZE]{};
      for(int i = 0; i < SIZE; i++)
        a[i] = m[i].count();
      double ave = 0;
      for(int i = 0; i < SIZE; i++)
        ave += a[i];
      ave /= SIZE;
      for(int i = 0; i < SIZE; i++)
        ret += (a[i] - ave) * (a[i] - ave);
      ret /= SIZE;
      return ret;
    }

    double H_var(){
      return this->h_var() + this->inv().trans().h_var();
    }

    double h_my(){
      double ret = 0;
      for(int i = 0; i < SIZE; i++){
        double tmp = 0;
        tmp += m[i].count();
        for(int j = 0; j < SIZE; j++)
          if(j != i)
            tmp += m[i][j];
        ret += tmp * tmp;
      }
      return ret;
    }

    double H_my(){
      return h_my() + this->inv().h_my();
    }


    double h_sq(){
      double ret = 0;
      for(int i = 0; i < SIZE; i++)
          ret += m[i].count() * m[i].count();
      return ret;
    }

    double h_sq_t(){
      double ret = 0;
        for(int i = 0; i < SIZE; i++){
          int cnt = 0;
          for(int j = 0; j < SIZE; j++)
            if(m[j][i] == 1)
              cnt++;
          ret += cnt * cnt;
        }
      return ret;
    }


    double H_sq(){
      return this->h_sq() + this->inv().h_sq_t();
    }

    double HT_sq(){
      return this->h_sq() + this->h_sq_t() + this->inv().h_sq()+ this->inv().h_sq_t();
    }


    double h_sum(){
        double ret = 0;
        for(int i = 0; i < SIZE; i++)
            ret += m[i].count();
        return ret;
    }
    

    double H_sum(){
        return this->h_sum() + this->inv().h_sum();
    }
    
    double h_prod(){
        double ret = 0;
        for(int i = 0; i < SIZE; i++)
            ret += log2(m[i].count());
        return ret;
    }

    double h_prod_t(){
        double ret = 0;
        for(int i = 0; i < SIZE; i++){
          int cnt = 0;
          for(int j = 0; j < SIZE; j++)
            if(m[j][i] == 1)
              cnt++;
          ret += log2(cnt);
        }
            
        return ret;
    }

    
    double H_prod(){
        return this->h_prod() + this->inv().h_prod();
    }
    
    double HT_prod(){
        return this->h_prod() + this->inv().h_prod() + this->h_prod_t() + this->inv().h_prod_t();
    }

    bool can_depthone(){
        int mx_w = 0;
        for(int i = 0; i < SIZE; i++)
            if(m[i].count() > 2)
                return false;
        int cnt_col[SIZE]{};
        for(int i = 0; i < SIZE; i++){
            int cnt = 0;
            for(int j = 0; j < SIZE; j++)
                if(m[j][i])
                    cnt++;
            if(cnt > 2)
                return false;
            cnt_col[i] = cnt;
        }

        int cnt_r[SIZE]{};
        int cnt_c[SIZE]{};
        for(int i = 0; i < SIZE; i++)
            for(int j = 0; j < SIZE; j++){
                if(m[i][j] && m[i].count() == 2 && cnt_col[j] == 2)
                    cnt_r[i]++, cnt_c[j]++;
                if(cnt_r[i] > 1 || cnt_c[j] > 1)
                    return false;
            }
        return true;
    }
    
};


typedef struct{
    int src;
    int tar;
    int type;//0: row; 1: col.
}Op;




void greedy0(Mat mat){

    Mat origin = mat;

    Mat inv = mat.inv();
    Mat tmp_mat;
    Mat tmp_inv;


    int depth = 0;
    double minm_cst = DBL_MAX; //large
    double tmp_cst;
    double minm_tie;
    double tmp_tie;
    
    

    vector<Op> select_list;
    vector<Op> r_layer, c_layer;
    vector<vector<Op> > r_layers, c_layers;
    bitset<SIZE> row_visi, col_visi;
    vector<Op> r_op, c_op;
    bool one = false;
    while(mat.cst_func(5) + inv.trans().cst_func(5) > costfun_min[6]){
        select_list.clear();
        minm_cst = max(mat.cst_func(5) + inv.trans().cst_func(5), mat.trans().cst_func(5) + inv.cst_func(5));
        // minm_cst = mat.cst_func(5) + inv.trans().cst_func(5);
        //row trans
        if(!one){

        for(int i = 0; i < SIZE; i++){
            if(row_visi[i] == 1)
                continue;
            for(int j = 0; j < SIZE; j++){
                if(row_visi[j] == 1 || j == i)
                    continue;
                tmp_mat = mat.row_i2j(i, j);
                tmp_inv = inv.col_i2j(j, i);
                tmp_cst = tmp_mat.cst_func(5) + tmp_inv.trans().cst_func(5);
                // tmp_tie = tmp_mat.trans().cst_func(11);

                if(tmp_cst < minm_cst + epsilon){
                    if(tmp_cst < minm_cst - epsilon){
                        select_list.clear();
                        Op op = {i, j, 0};
                        select_list.push_back(op);
                        minm_cst = tmp_cst;
                        // minm_tie = tmp_tie;
                    }
                    // else if(tmp_tie < minm_tie){
                    //     select_list.clear();
                    //     Op op = {i, j, 0};
                    //     select_list.push_back(op);
                    //     minm_tie = tmp_tie;
                    // }
                    else{
                        Op op = {i, j, 0};
                        select_list.push_back(op);                      
                    }
                }
            }
        }

        }
        
        //col trans
        for(int i = 0; i < SIZE; i++){
            if(col_visi[i] == 1)
                continue;
            for(int j = 0; j < SIZE; j++){
                if(col_visi[j] == 1 || j == i)
                    continue;
                tmp_mat = mat.col_i2j(i, j);
                tmp_inv = inv.row_i2j(j, i);
                tmp_cst = tmp_mat.trans().cst_func(5) + tmp_inv.cst_func(5);
                // tmp_tie = tmp_mat.cst_func(11);//trans与否？
                if(tmp_cst < minm_cst + epsilon){
                    if(tmp_cst < minm_cst - epsilon){
                        select_list.clear();
                        Op op = {i, j, 1};
                        select_list.push_back(op);
                        minm_cst = tmp_cst;
                        // minm_tie = tmp_tie;
                    }
                    // else if(tmp_tie < minm_tie){
                    //     select_list.clear();
                    //     Op op = {i, j, 1};
                    //     select_list.push_back(op);
                    //     minm_tie = tmp_tie;
                    // }
                    else{
                        Op op = {i, j, 1};
                        select_list.push_back(op);                      
                    }
                }
            }
        }
        
        if(select_list.size() == 0){
        //   cout << RANK << " depth " << depth << endl;
        //   mat.print();
            if(r_layer.size()){
                r_layers.push_back(r_layer);
                // r_op.insert(r_op.end(), r_layer.begin(), r_layer.end());
                // cout << RANK << r_layer.size() << ' ';
                r_layer.clear();
                row_visi = 0;
            }
            if(c_layer.size()){
                c_layers.push_back(c_layer);
                // cout << RANK << c_layer.size() << ' ';
                // c_op.insert(c_op.end(), c_layer.begin(), c_layer.end());
                c_layer.clear();
                col_visi = 0;
            }
            // cout << endl;
            if(mat.can_depthone()){
                one = true;
                // mat.print();
                // int cnt_col[SIZE]{};
                // for(int i = 0; i < SIZE; i++)
                //     for(int j = 0; j < SIZE; j++)
                //         if(mat.m[i][j])
                //             cnt_col[j]++;
                // for(int i = 0; i < SIZE; i++)
                //     for(int j = 0; j < SIZE; j++){
                //         if(mat.m[i].count() == 2 && cnt_col[j] == 2){
                //             for(int k = 0; k < SIZE; k++)
                //                 if(k != j && mat.m[i][k]){
                //                     mat = mat.col_i2j(k, j);
                //                     c_layer.push_back({k, j, 1});
                //                     c_op.push_back({k, j, 1});
                //                 }
                //         }
                //     }
                // mat.print();
                
                // depth++;
                // break;
            }
        }
        else{
            int rd = rand() % select_list.size();
            auto [c, t, t01] = select_list[rd];
            if(t01 == 0){
                mat = mat.row_i2j(c, t);
                inv = inv.col_i2j(t, c);
                r_layer.push_back({c, t, 0});
                r_op.push_back({c, t, 0});
                if(row_visi == 0)
                    depth++;
                row_visi[c] = 1;
                row_visi[t] = 1;
            }
            else{
                mat = mat.col_i2j(c, t);
                inv = inv.row_i2j(t, c);
                c_layer.push_back({c, t, 1});
                c_op.push_back({c, t, 1});
                if(col_visi == 0)
                    depth++;
                col_visi[c] = 1;
                col_visi[t] = 1;
            }
            // minm_cst = mat.cst_func(5) + inv.trans().cst_func(5);
        }
        
        if(depth > 100){
            cout << "d too large" << endl;
            mat.print();
            break;
        }
    }



    if(r_layer.size()){
        // r_op.insert(r_op.end(), r_layer.begin(), r_layer.end());
        // cout << RANK << r_layer.size() << ' ';
        r_layers.push_back(r_layer);
        r_layer.clear();
        row_visi = 0;
    }



    if(c_layer.size()){
        // cout << RANK << c_layer.size() << ' ';
        c_layers.push_back(c_layer);
        // c_op.insert(c_op.end(), c_layer.begin(), c_layer.end());
        c_layer.clear();
        col_visi = 0;
    }
    // cout << endl;
    // cout << RANK << "depth is " << depth << endl;
    // mat.print();
    
    
    // 打印reduce结果
    // cout << "aes reduce" << endl;
    Mat reduce = origin;
    int sz = 0;//cnot count
    for(auto [c, t, _]: r_op)
        reduce = reduce.row_i2j(c, t), sz++;
    for(auto [c, t, _]: c_op)
        reduce = reduce.col_i2j(c, t), sz++;
    // cout << RANK << "  : " << endl;
    // reduce.print();
    // cout << endl;

    bool ok = true;
    for(int i = 0; i < SIZE; i++)
        if(reduce.m[i] != mat.m[i]){
            ok = false;
            cout << i << " th row diff" << endl;
        }
    if(depth > minm || (depth == minm && sz >= sz_minm) || !ok)
        return;
    
    minm = depth;
    sz_minm = sz;
    cout << CHOICE << " " << RANK << "  d:" << minm << "  sz:" << sz_minm << " " << FILENAME << " OPT:" << OPT << endl;

    int per[SIZE];
    for(int i = 0; i < SIZE; i++)
        for(int j = 0; j < SIZE; j++)
            if(mat.m[i][j] == 1)
                per[i] = j;

    vector<Op> seq;
    vector<vector<Op> > layers;

    for(int i = 0; i < c_op.size(); i++){
        auto [t, c, _] = c_op[i];
        seq.push_back({c, t, 1});
    }
    for(int i = r_op.size() - 1; i >= 0 ; i--){
        auto [c, t, _] = r_op[i];
        seq.push_back({per[c], per[t], 0});
    }
    


    // cout << FILENAME << RANK << "c_layers.size()" << c_layers.size() << "r_layers.size()" << r_layers.size() << endl;
    for(auto& l: c_layers){
        vector<Op> nl;
        for(auto [t, c, _]: l){
            nl.push_back({c, t, 1});
        }
        layers.push_back(nl);
    }
    reverse(r_layers.begin(), r_layers.end());//列1-d， 行d到1
    for(auto& l: r_layers){
        vector<Op> nl;
        for(auto [c, t, _]: l){
            nl.push_back({per[c], per[t], 0});
        }
        layers.push_back(nl);
    }
    string mat_str = M_string;

    

    



    

    string res = "res/";
    ofstream f_layer(res + "layers_" + RANK + string(FILENAME) , ios::app);
    for(auto& l: layers){
        for(auto [t, c, _]: l){
            f_layer << t << ' ' << c << ' ';
        }
        f_layer << endl;
    }
    f_layer << "CNOT: " << seq.size() << " depth: " << layers.size() << endl << endl;
    f_layer.close();


    
    ofstream fout(res + "seq_" + RANK + FILENAME , ios::app);
    for(int i = 0; i < seq.size(); i++)
        fout << seq[i].src << ' ' << seq[i].tar << ' ' << seq[i].type << endl;
    fout << "CNOT: " << seq.size() << endl << endl;


    // ofstream fnote("note.txt", ios::app);
    // fnote << RANK << FILENAME << " " << layers.size() << " depth " << seq.size() << endl;
    // fnote.close(); 


    Mat verify;
    for(int i = 0; i < SIZE; i++)
        verify.m[i][i] = 1;
    
    for(int i = 0; i < seq.size(); i++)
        verify = verify.row_i2j(seq[i].src, seq[i].tar);
    for(int i = 0; i < SIZE; i++){
        bool ok = false;
        for(int j = 0; j < SIZE; j++)
            if(verify.m[i] == origin.m[j])
                ok = true;
        if(!ok)
            cout << "row " << i << " does not exist." << endl;
    }

    // verify.print();

    
    fout.close();

}


void greedy_origin(Mat mat){

    Mat origin = mat;

    Mat inv = mat.inv();
    Mat tmp_mat;
    Mat tmp_inv;


    int depth = 0;
    double minm_cst = DBL_MAX; //large
    double tmp_cst;
    double minm_tie;
    double tmp_tie;
    
    

    vector<Op> select_list;
    vector<Op> r_layer, c_layer;
    vector<vector<Op> > r_layers, c_layers;
    bitset<SIZE> row_visi, col_visi;
    vector<Op> r_op, c_op;
    bool one = false;
    while(mat.cst_func(3) + inv.cst_func(3) > costfun_min[4]){
        select_list.clear();
        minm_cst = mat.cst_func(3) + inv.cst_func(3);
        // minm_cst = mat.cst_func(5) + inv.trans().cst_func(5);
        //row trans
        if(!one){

        for(int i = 0; i < SIZE; i++){
            if(row_visi[i] == 1)
                continue;
            for(int j = 0; j < SIZE; j++){
                if(row_visi[j] == 1 || j == i)
                    continue;
                tmp_mat = mat.row_i2j(i, j);
                tmp_inv = inv.col_i2j(j, i);
                tmp_cst = tmp_mat.cst_func(3) + tmp_inv.cst_func(3);
                // tmp_tie = tmp_mat.trans().cst_func(11);

                if(tmp_cst < minm_cst + epsilon){
                    if(tmp_cst < minm_cst - epsilon){
                        select_list.clear();
                        Op op = {i, j, 0};
                        select_list.push_back(op);
                        minm_cst = tmp_cst;
                        // minm_tie = tmp_tie;
                    }
                    // else if(tmp_tie < minm_tie){
                    //     select_list.clear();
                    //     Op op = {i, j, 0};
                    //     select_list.push_back(op);
                    //     minm_tie = tmp_tie;
                    // }
                    else{
                        Op op = {i, j, 0};
                        select_list.push_back(op);                      
                    }
                }
            }
        }

        }
        
        //col trans
        for(int i = 0; i < SIZE; i++){
            if(col_visi[i] == 1)
                continue;
            for(int j = 0; j < SIZE; j++){
                if(col_visi[j] == 1 || j == i)
                    continue;
                tmp_mat = mat.col_i2j(i, j);
                tmp_inv = inv.row_i2j(j, i);
                tmp_cst = tmp_mat.cst_func(3) + tmp_inv.cst_func(3);
                // tmp_tie = tmp_mat.cst_func(11);//trans与否？
                if(tmp_cst < minm_cst + epsilon){
                    if(tmp_cst < minm_cst - epsilon){
                        select_list.clear();
                        Op op = {i, j, 1};
                        select_list.push_back(op);
                        minm_cst = tmp_cst;
                        // minm_tie = tmp_tie;
                    }
                    // else if(tmp_tie < minm_tie){
                    //     select_list.clear();
                    //     Op op = {i, j, 1};
                    //     select_list.push_back(op);
                    //     minm_tie = tmp_tie;
                    // }
                    else{
                        Op op = {i, j, 1};
                        select_list.push_back(op);                      
                    }
                }
            }
        }
        
        if(select_list.size() == 0){
        //   cout << RANK << " depth " << depth << endl;
        //   mat.print();
            if(r_layer.size()){
                r_layers.push_back(r_layer);
                // r_op.insert(r_op.end(), r_layer.begin(), r_layer.end());
                // cout << RANK << r_layer.size() << ' ';
                r_layer.clear();
                row_visi = 0;
            }
            if(c_layer.size()){
                c_layers.push_back(c_layer);
                // cout << RANK << c_layer.size() << ' ';
                // c_op.insert(c_op.end(), c_layer.begin(), c_layer.end());
                c_layer.clear();
                col_visi = 0;
            }
            // cout << endl;
            // if(mat.can_depthone()){
            //     one = true;
            //     // mat.print();
            //     // int cnt_col[SIZE]{};
            //     // for(int i = 0; i < SIZE; i++)
            //     //     for(int j = 0; j < SIZE; j++)
            //     //         if(mat.m[i][j])
            //     //             cnt_col[j]++;
            //     // for(int i = 0; i < SIZE; i++)
            //     //     for(int j = 0; j < SIZE; j++){
            //     //         if(mat.m[i].count() == 2 && cnt_col[j] == 2){
            //     //             for(int k = 0; k < SIZE; k++)
            //     //                 if(k != j && mat.m[i][k]){
            //     //                     mat = mat.col_i2j(k, j);
            //     //                     c_layer.push_back({k, j, 1});
            //     //                     c_op.push_back({k, j, 1});
            //     //                 }
            //     //         }
            //     //     }
            //     // mat.print();
                
            //     // depth++;
            //     // break;
            // }
        }
        else{
            int rd = rand() % select_list.size();
            auto [c, t, t01] = select_list[rd];
            if(t01 == 0){
                mat = mat.row_i2j(c, t);
                inv = inv.col_i2j(t, c);
                r_layer.push_back({c, t, 0});
                r_op.push_back({c, t, 0});
                if(row_visi == 0)
                    depth++;
                row_visi[c] = 1;
                row_visi[t] = 1;
            }
            else{
                mat = mat.col_i2j(c, t);
                inv = inv.row_i2j(t, c);
                c_layer.push_back({c, t, 1});
                c_op.push_back({c, t, 1});
                if(col_visi == 0)
                    depth++;
                col_visi[c] = 1;
                col_visi[t] = 1;
            }
            // minm_cst = mat.cst_func(5) + inv.trans().cst_func(5);
        }
        
        if(depth > 100){
            cout << "d too large" << endl;
            break;
        }
    }



    if(r_layer.size()){
        // r_op.insert(r_op.end(), r_layer.begin(), r_layer.end());
        // cout << RANK << r_layer.size() << ' ';
        r_layers.push_back(r_layer);
        r_layer.clear();
        row_visi = 0;
    }



    if(c_layer.size()){
        // cout << RANK << c_layer.size() << ' ';
        c_layers.push_back(c_layer);
        // c_op.insert(c_op.end(), c_layer.begin(), c_layer.end());
        c_layer.clear();
        col_visi = 0;
    }
    // cout << endl;
    // cout << RANK << "depth is " << depth << endl;
    // mat.print();
    Mat reduce = origin;
    int sz = 0;//cnot count
    for(auto [c, t, _]: r_op)
        reduce = reduce.row_i2j(c, t), sz++;
    for(auto [c, t, _]: c_op)
        reduce = reduce.col_i2j(c, t), sz++;
    // cout << RANK << "  : " << endl;
    // reduce.print();
    // cout << endl;

    bool ok = true;
    for(int i = 0; i < SIZE; i++)
        if(reduce.m[i] != mat.m[i]){
            ok = false;
            cout << i << " th row diff" << endl;
        }
    if(depth > minm || (depth == minm && sz >= sz_minm) || !ok)
        return;
    
    minm = depth;
    sz_minm = sz;
    cout << CHOICE << " " << RANK << "  d:" << minm << "  sz:" << sz_minm << " " << FILENAME << " OPT:" << OPT << endl;

    int per[SIZE];
    for(int i = 0; i < SIZE; i++)
        for(int j = 0; j < SIZE; j++)
            if(mat.m[i][j] == 1)
                per[i] = j;

    vector<Op> seq;
    vector<vector<Op> > layers;

    for(int i = 0; i < c_op.size(); i++){
        auto [t, c, _] = c_op[i];
        seq.push_back({c, t, 1});
    }
    for(int i = r_op.size() - 1; i >= 0 ; i--){
        auto [c, t, _] = r_op[i];
        seq.push_back({per[c], per[t], 0});
    }
    


    // cout << FILENAME << RANK << "c_layers.size()" << c_layers.size() << "r_layers.size()" << r_layers.size() << endl;
    for(auto& l: c_layers){
        vector<Op> nl;
        for(auto [t, c, _]: l){
            nl.push_back({c, t, 1});
        }
        layers.push_back(nl);
    }
    reverse(r_layers.begin(), r_layers.end());//列1-d， 行d到1
    for(auto& l: r_layers){
        vector<Op> nl;
        for(auto [c, t, _]: l){
            nl.push_back({per[c], per[t], 0});
        }
        layers.push_back(nl);
    }
    string mat_str = M_string;

    

    



    

    string res = "res/";
    ofstream f_layer(res + "layers_" + RANK + string(FILENAME) , ios::app);
    for(auto& l: layers){
        for(auto [t, c, _]: l){
            f_layer << t << ' ' << c << ' ';
        }
        f_layer << endl;
    }
    f_layer << "CNOT: " << seq.size() << " depth: " << layers.size() << endl << endl;
    f_layer.close();


    
    ofstream fout(res + "seq_" + RANK + FILENAME , ios::app);
    for(int i = 0; i < seq.size(); i++)
        fout << seq[i].src << ' ' << seq[i].tar << ' ' << seq[i].type << endl;
    fout << "CNOT: " << seq.size() << endl << endl;


    // ofstream fnote("note.txt", ios::app);
    // fnote << RANK << FILENAME << " " << layers.size() << " depth " << seq.size() << endl;
    // fnote.close(); 


    Mat verify;
    for(int i = 0; i < SIZE; i++)
        verify.m[i][i] = 1;
    
    for(int i = 0; i < seq.size(); i++)
        verify = verify.row_i2j(seq[i].src, seq[i].tar);
    for(int i = 0; i < SIZE; i++){
        bool ok = false;
        for(int j = 0; j < SIZE; j++)
            if(verify.m[i] == origin.m[j])
                ok = true;
        if(!ok)
            cout << "row " << i << " does not exist." << endl;
    }

    // verify.print();

    
    fout.close();
}


Mat read_mat_xiang(){
    string folder1 = "Implementations/";
    string s_file = folder1 + FILENAME;
    fstream fin(s_file);
    string s;
    vector<Op> SEQ;
    Mat mat;
    getline(fin, s);
    for(int i = 0; i < SIZE; i++){
        getline(fin, s);
        for(int j = 0; j < SIZE; j++)
            mat.m[i][j] = s[j] - '0';
    }



    while(getline(fin, s)){
        if(s.size() >= 1 && s[0] == 'X')
            break;
    }

    while(getline(fin, s)){
        vector<int> nums;
        int x = 0;
        bool isnum = 0;
        for(char c: s){
            if(c >= '0' && c <= '9'){
                isnum = true;
                x = x * 10 + c - '0';
            }
            else{
                if(isnum)
                    nums.push_back(x);
                isnum = false;
                x = 0;
            }
        }
        if(isnum)
            nums.push_back(x);
        bool y = -1;
        if(nums.size() == 7){//不看flag
            y = nums[3];
        }
        SEQ.push_back({nums[2], nums[0], y});
    }
    fin.close();
    return mat;
}




Mat read_mat(){
  int x, y, z, tmp;
  string folder = "matrices_bp_format/";
  fstream fin(folder + FILENAME);
  fin >> x >> y >> z;
  Mat mat;
  for(int i = 0; i < SIZE; i++)
    for(int j = 0; j < SIZE; j++){
      fin >> tmp;
      mat.m[i][j] = tmp;
    }
  return mat;
}

void run_mpi(int argc, char** argv){
    MPI_Init(&argc, &argv);
    int world_size;
    MPI_Comm_size(MPI_COMM_WORLD, &world_size);
    int world_rank;
    MPI_Comm_rank(MPI_COMM_WORLD, &world_rank);   
    int a;
    srand(time(NULL) + (unsigned long)(&a)); 
    
    
    RANK = to_string(world_rank) + "____";
    string mx_duration(argv[1]);
    int mx_dr = stoi(mx_duration);
    clock_t st, en;
    st = clock();
    // fstream fin("aes.txt");
    // string s;
    // Mat aes;
    // int i = 0;
    // while(getline(fin, s)){
    //     for(int j = 0; j < SIZE; j++)
    //         aes.m[i][j] = s[j] - '0';
    //     i++;
    // }
    // fin.close();


    // aes.print();
    // Mat aesinv = aes.inv();

    // greedy0_tie(aes);
    // Mat mat = read_mat_xiang();
    string res = "res/" + string(FILENAME);
    fstream fopt(res + string(FILENAME));
    int tmp;
    if((fopt >> tmp)){
        minm = tmp;
    }
    fopt.close();

    Mat mat = read_mat();
    long long cnt = 0;
    while(true){
        cnt++;
        if(cnt % 100 == 0)
            cout << RANK << "cnt " << cnt << endl;
        greedy0(mat);
        en = clock();
        if((double)(en - st) / CLOCKS_PER_SEC > mx_dr)
            break;
    }
    // ofstream fnotes("notes.txt", ios::app);
    // fnotes << CHOICE << " " << RANK << "  " << minm << "  " << FILENAME << " OPT:" << OPT << "count " << cnt << endl;
    // fnotes.close();

    cout << CHOICE << " " << RANK << "  d:" << minm << "  sz:" << sz_minm << " " << FILENAME << " OPT:" << OPT << "count " << cnt << endl;

	string record = "record/" + RANK + string(FILENAME);
	ofstream f(record, ios::app);
    f << minm << "  " << sz_minm << endl;
    f.close();

    MPI_Finalize();
}






void blossom_greedy(Mat mat){
    Mat origin = mat;
    Mat tmp_mat;
    Mat inv_mat;
    Mat inv_tmp_mat;
    int depth = 0;
    double minm_cst = DBL_MAX; //large
    double tmp_cst;
    
    

    vector<Op> select_list;
    vector<Op> r_layer, c_layer;
    vector<vector<Op> > r_layers, c_layers;
    bitset<SIZE> row_visi, col_visi;
    vector<Op> r_op, c_op;
    while(mat.cst_func(FUNC) > costfun_min[FUNC]){
        minm_cst = mat.cst_func(FUNC);
        inv_mat = mat.inv();
        
        r_layer.clear();
        ma.init(SIZE);
        //row trans
        for(int i = 0; i < SIZE; i++){
            for(int j = 0; j < SIZE; j++){
                if(j == i)
                    continue;
                
                tmp_mat = mat.row_i2j(i, j);
                inv_tmp_mat = inv_mat.col_i2j(j, i);
                int reduce_ij = mat.cst_func(1) + inv_mat.cst_func(1) - (tmp_mat.cst_func(1) + inv_tmp_mat.cst_func(1));//to opt

                tmp_mat = mat.row_i2j(j, i);
                inv_tmp_mat = inv_mat.col_i2j(i, j);
                int reduce_ji = mat.cst_func(1) + inv_mat.cst_func(1) - (tmp_mat.cst_func(1) + inv_tmp_mat.cst_func(1));

                int mx = max(reduce_ij, reduce_ji);
                if(mx > 0)
                    ma.add_edge(i, j, mx);
                
            }
        }

        int w1 = ma.solve().first;

        bitset<SIZE> visited;
        for(int i = 1; i <= SIZE; i++){
            if(ma.match[i] > 0 && visited[i] == 0){
                int j = ma.match[i];
                visited[i] = 1, visited[j] = 1;

                tmp_mat = mat.row_i2j(i, j);
                inv_tmp_mat = inv_mat.col_i2j(j, i);
                int reduce_ij = mat.cst_func(1) + inv_mat.cst_func(1) - (tmp_mat.cst_func(1) + inv_tmp_mat.cst_func(1));

                tmp_mat = mat.row_i2j(j, i);
                inv_tmp_mat = inv_mat.col_i2j(i, j);
                int reduce_ji = mat.cst_func(1) + inv_mat.cst_func(1) - (tmp_mat.cst_func(1) + inv_tmp_mat.cst_func(1));
                if(reduce_ij > reduce_ji){
                    mat = mat.row_i2j(i, j);
                    inv_mat = inv_mat.col_i2j(j, i);
                    r_layer.push_back({i, j, 0});
                }
                else if(reduce_ij <= reduce_ji){//equal, random?
                    mat = mat.row_i2j(j, i);
                    inv_mat = inv_mat.col_i2j(i, j);
                    r_layer.push_back({j, i, 0});
                }
            }
        }
        r_layers.push_back(r_layer);
        depth++;

        cout << depth << " depth with " << mat.cst_func(1) << " cst func" << endl;
        for(auto [c, t, _]: r_layer){
            cout << c << "  " << t << " " << endl;
        }
        cout << endl;

        
        
        if(depth > 60){
            cout << "d too large" << endl;
            break;
        }
    }
    
    cout << endl;
    cout << "depth is " << depth << endl;
    mat.print();


    cout << "aes reduce" << endl;
    Mat reduce = origin;
    for(auto [c, t, _]: r_op)
        reduce = reduce.row_i2j(c, t);
    for(auto [c, t, _]: c_op)
        reduce = reduce.col_i2j(c, t);
    reduce.print();
    for(int i = 0; i < SIZE; i++)
        if(reduce.m[i] != mat.m[i])
            cout << i << " th row diff" << endl;
    
    int per[SIZE];
    for(int i = 0; i < SIZE; i++)
        for(int j = 0; j < SIZE; j++)
            if(mat.m[i][j] == 1)
                per[i] = j;

    vector<Op> seq;
    vector<vector<Op> > layers;
    for(auto& l: c_layers){
        vector<Op> nl;
        for(auto [t, c, _]: l){
            nl.push_back({c, t, 1});
        }
        layers.push_back(nl);
    }
    for(auto& l: r_layers){
        vector<Op> nl;
        for(auto [c, t, _]: l){
            nl.push_back({per[c], per[t], 0});
        }
        layers.push_back(nl);
    }
    string mat_str = M_string;
    ofstream f_layer(mat_str + "layers.txt", ios::app);
    for(auto& l: layers){
        for(auto [t, c, _]: l){
            f_layer << t << ' ' << c << ' ';
        }
        f_layer << endl;
    }
    f_layer << endl;
    f_layer.close();



    for(int i = 0; i < c_op.size(); i++){
        auto [t, c, _] = c_op[i];
        seq.push_back({c, t, 1});
    }
    for(int i = r_op.size() - 1; i >= 0 ; i--){
        auto [c, t, _] = r_op[i];
        seq.push_back({per[c], per[t], 0});
    }
        
    ofstream fout(mat_str + "seq.txt", ios::app);
    for(int i = 0; i < seq.size(); i++)
        fout << seq[i].src << ' ' << seq[i].tar << ' ' << seq[i].type << endl;
    fout << endl;

    Mat verify;
    for(int i = 0; i < SIZE; i++)
        verify.m[i][i] = 1;
    
    for(int i = 0; i < seq.size(); i++)
        verify = verify.row_i2j(seq[i].src, seq[i].tar);
    for(int i = 0; i < SIZE; i++){
        bool ok = false;
        for(int j = 0; j < SIZE; j++)
            if(verify.m[i] == origin.m[j])
                ok = true;
        if(!ok)
            cout << "row " << i << " does not exist." << endl;
    }

    verify.print();

    fout.close();

}


int main(int argc, char** argv)
{
    // bitset<SIZE> a("00000000000000000000000000000001");
    // cout << a << endl;
    // a = 0;
    // // a[0] = 1;
    // cout << a << endl;
    // bitset<SIZE> m[SIZE];
    // string s = "00000000000000000000000000000001";
    // m[0] = bitset<SIZE>(s);
    // cout << m[0] << endl;
    // Mat a;
    // a.m[0][0] = 0;
    // a.m[0][1] = 1;
    // a.m[1][0] = 1;
    
    // cout << a.h_sum() << ' ' << a.H_sum() << ' ' << a.h_prod() << ' ' << a.H_prod() << endl;
    // // cout << "      " << (a.m[0][0] == 1) << endl;
    // bitset<SIZE> b[SIZE];
    // b[0][0] = 1;
    // b[0][1] = 1;
    // b[1][1] = 1;
    
    // a.print();
    // Mat c = a.inv();
    // c.print();    
    
    srand(time(NULL));
    


    run_mpi(argc, argv);
    // blossom_greedy(aes);
 
    
    
    return 0;
}
